{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will show how to calculate differences between a DEM and ICESat2 points using the `reference_dem` class.\n",
    "\n",
    "First, import `coregistration` and create the DEM object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import the module via an explicit path directly for now\n",
    "import importlib.util\n",
    "spec = importlib.util.spec_from_file_location(\"coregistration\", \"/home/jovyan/Assimilation/scripts/coregistration.py\")\n",
    "cor = importlib.util.module_from_spec(spec)\n",
    "spec.loader.exec_module(cor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_file_path = '/home/jovyan/data/reference_dem_clip.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = cor.reference_dem(dem_file_path)\n",
    "ds.calculate_bounding_box(4326)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, create a GeoDataFrame of ATL06 points. This will be done by our scripts later but is provided explicitly as an example for now. \n",
    "\n",
    "First download an example file using the scripts provided by Tian Li.\n",
    "\n",
    "Then, create a GeoDataFrame from one of the files (e.g. the `gt2r` beam points on 2018-10-19):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "from pyproj import CRS\n",
    "\n",
    "example_file = '/home/jovyan/data/ATL06_Processed/processed_ATL06_20181019231110_03260102_003_01_gt2r.h5'\n",
    "vnames = ['x', 'y', 'h_li']\n",
    "i2_data = atl06_lib.read_h5(example_file, vnames)\n",
    "\n",
    "points = []            \n",
    "for p in range(np.shape(i2_data)[0]):\n",
    "    points.append({\n",
    "        'x': i2_data[p,0],\n",
    "        'y': i2_data[p,1],\n",
    "        'h': i2_data[p,2]})\n",
    "\n",
    "df = pd.DataFrame.from_dict(points)\n",
    "\n",
    "geometry = [Point(xyz) for xyz in zip(df['x'], df['y'], df['h'])]        \n",
    "gdf = gpd.GeoDataFrame(df, geometry=geometry, crs=CRS('EPSG:'+str(ds.epsg)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the `colocate_icesat2_dem_points(icesat2_gdf)` can be called to find the difference between the provided GeoDataFrame of icesat2 points and the DEM points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 x             y     dem_elev     is2_elev    dem-is2\n",
      "1    586319.497387  5.395103e+06  1039.610107  1038.600952   1.009155\n",
      "2    586317.188128  5.395124e+06  1054.310059  1039.262085  15.047974\n",
      "3    586314.805116  5.395145e+06  1040.072998  1040.090332  -0.017334\n",
      "4    586312.578000  5.395165e+06  1042.106812  1040.767700   1.339111\n",
      "5    586310.466515  5.395184e+06  1041.342896  1041.149780   0.193115\n",
      "..             ...           ...          ...          ...        ...\n",
      "139  585518.317682  5.402210e+06  2480.963867  2485.448730  -4.484863\n",
      "140  585516.129457  5.402230e+06  2482.430664  2486.523926  -4.093262\n",
      "141  585513.934961  5.402250e+06  2484.899658  2488.912842  -4.013184\n",
      "142  585486.745317  5.402489e+06  2606.960938  2607.850586  -0.889648\n",
      "143  585484.533362  5.402509e+06  2611.922852  2613.419678  -1.496826\n",
      "\n",
      "[143 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "ds.colocate_icesat2_dem_points(gdf)\n",
    "print(ds.dem_is2_comparison)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the comparison of points are stored in the `dem_reference` object"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
